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The solution-space structure of the 3-Satisfiability Problem (3-SAT) is studied as a function of the control 
parameter a (ratio of number of clauses to the number of variables) using numerical simulations. For this 
purpose, one has to sample the solution space with uniform weight. It is shown here that standard stochastic 
local-search (SLS) algorithms like "ASAT" and "MCMCMC" (also known as "parallel tempering") exhibit a 
sampling bias. Nevertheless, unbiased samples of solutions can be obtained using the "ballistic-networking 
approach", which is introduced here. It is a generalization of "ballistic search" methods and yields also a cluster 
structure of the solution space. 

As application, solutions of 3-SAT instances are generated using ASAT plus ballistic networking. The nu- 
merical results are compatible with a previous analytic prediction of a simple solution-space structure for small 
values of a and a transition to a clustered phase at oc c ~ 3.86, where the solution space breaks up into several 
non-negligible clusters. Furthermore, in the thermodynamic limit there are, for values of a close to the SAT- 
UNSAT transition oc s ~ 4.267, always clusters without any frozen variables. This may explain why some SLS 
algorithms are able to solve very large 3-SAT instances close to the SAT-UNSAT transition. 
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I. INTRODUCTION 

The application of notions, analytical approaches and nu- 
merical algorithms from statistical mechanics has lead to a 
better understanding [1-3] of NP-hard optimization problems 
HHH]. One main underlying question is, why these optimiza- 
tion problems are computationally hard. This means no fast 
algorithms are available, where the running times increase 
only polynomially with the problem size. The progress of 
gaining insight into this phenomenon has been considerable in 
particular for the typical-case complexity, where ensembles of 
random instances are studied as a function of control param- 
eters. These ensembles often exhibit phase transitions where 
changes of the effective "hardness" of the problem can be ob- 
served. Often, these transitions are connected to changes of 
the structure of the solution space, comparable to energy land- 
scapes in physics. In particular, one is interested in the ques- 
tion, how the change of the solution-space structures has an in- 
fluence on the performance of exact and stochastic algorithms. 
For example, for the vertex-cover problem, which is defined 
on graphs, a clustering transition has been found analytically 
|0] and numerically J3H] when increasing the edge density of 
Erdos-Renyi random graphs. This transition coincides with a 
change of the typical-case complexity from polynomial to ex- 
ponential |9||. For other optimization problems, the situation 
is less clear, as for the satisfiability problem (SAT), which we 
study in this work. 

As we will explain, exact enumeration of solutions works 
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well in one region of the phase diagram, close to the SAT- 
UNSAT phase transition (see below), whereas Monte Carlo 
approaches perform well in the opposite part of the phase dia- 
gram, away from the SAT-UNSAT transitions. Unfortunately, 
the clustering transition is located right between these extreme 
parts, hence numerically difficult to study. We use a stochas- 
tic algorithm in combination with a correction of the sampling 
bias introduced by the stochastic algorithm to study the clus- 
tering phenomena. 

The outline of the paper is as follows. In the second section, 
we give the necessary background on SAT and on clustering 
of solution landscapes. In the third section, we briefly explain 
the algorithms we use to sample solutions and show that they 
exhibit a bias. Next, we introduce ballistic networking and re- 
lated methods, which we use to correct for the bias. In section 
five, we show the results we have obtained for random 3-SAT. 
Finally, we provide a conclusion and an outlook. 



H. BACKGROUND 

A. Satisfiability 

Satisfiability is one of the fundamental problems of com- 
puter science, and has attracted a lot of attention over the 
past years, also by physicists, due to its similarity to spin- 
glass problems. It is the first problem proven to belong to the 
class of NP-complete problems HlOfl . a class of problems for 
which no algorithm has been found yet that exhibits a poly- 
nomial worst-case running time as a function of the problem 
size. Therefore it is still a challenge to find algorithms which 
perform well on typical instances and to understand the un- 
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derlying structure of the solution space which may hinder the 
performance of algorithms. 

Satisfiability belongs to the class of constraint satisfaction 
problems Ultl : Given N Boolean variables x, =0,1 and a 
Boolean formula F describing a set of constraints, each of 
which forbids a certain assignment of values to some of the 
variables, you are to decide whether F can be satisfied, i. e., 
whether there is an assignment x = {x\,... ,jc#) such that all 
constraints are fulfilled simultaneously. In the K-SAT formu- 
lation, F is given in conjunctive normal form, 

M 

F= A (CV/fV.-.V/f), 

711=1 

which describes a logical conjunction of M constraints 
(clauses) C m each containing a disjunction of K literals I™ 
(m = l,...,M;k= 1, . . . ,K) which are either a variable x, or a 
logically negated variable ~x\. 

A certain assignment of values to the variables is called 
a configuration in the following. If a configuration satisfies 
all clauses in F it is called a solution. In the random K- 
SAT ensemble each clause is chosen randomly and uniformly 
amongst the 2 K (^) possible combinations in which no vari- 
able appears twice. 

One defines a control parameter a = M/N which is the 
number of clauses M divided by the number of variables N. 
For low a the problem is typically satisfiable whereas for high 
values of a there typically is no solution lfl2[[l3ll . It has been 
proven rigorously IU4I1 that the transition between the satisfi- 
able and the unsatisfiable phase becomes sharp for ,/V — > °°. 
Whilst the position of the threshold for K = 2 is known ex- 
actly 11511 . for larger K there are only numerical estimates. In 
this paper we will stick to the K = 3 case, where every clause 
contains exactly three literals. The satisfiability transition is 
located in this case at a s = 4.267 lfl6ll . 



B. Cluster phenomena 

In addition to the SAT-UNSAT transition, analytical cal- 
culations lfl7i fl8ll give rise to evidence that there are further 
("structural") phase transitions which refer to the formation of 
disconnected clusters of solutions for high values of the con- 
trol parameter a in the satisfiable phase. Formally, clusters 
in constraint satisfaction problems can be defined as extremal 
Gibbs measures which gives the following picture for Satis- 
fiability: For small values of a all solutions are contained in 
one connected component (cluster). When a grows, more and 
more solutions disappear so that at some point the cluster de- 
composes into smaller clusters which initially, up to a thresh- 
old a c i, make up only an exponentially small fraction of all 
solutions, whereas above a c i many clusters contribute to the 
statisical behavior. Above a higher critical value oc c we en- 
ter another type of clustered phase which is dominated by a 
small number of large clusters. The case of 3-SAT is special, 
as here a ( / = a 6 , i. e., we directly enter the phase dominated 
by few clusters. The position of the dynamical threshold to 
the clustered phase is predicted to be at oc c « 3.86 11711 . 



This value is compatible with recent numerical results lfl9tl . 
where the cluster structure was investigated using the detec- 
tion of community structures. Unfortunately, the sampling 
was performed using an algorithm, which does not exhibit 
uniform sampling of the solutions, see below. Anyway, there 
is no general rule how to translate the formal definition of 
clusters, which holds in the thermodynamic limit, to finite sys- 
tem sizes, hence other approaches than community structures 
are possible. For numerical studies often a very appealing 
approach is used, where a cluster is defined as the connected 
components in a graph where each solution is represented by a 
vertex and edges connect solutions differing in only one vari- 
able. This definition of a cluster will be used in this work as 
well. For every two solutions belonging to the same cluster 
there is therefore a "path" of configurations which all solve 
the SAT instance at hand. Unfortunately, this path can be long 
and peppered with many dead ends or loops which makes it 
very difficult to decide whether two configurations belong to 
the same cluster. The main problem when discussing clusters 
in high-dimensional discrete solution spaces like that of Satis- 
fiability is that one is tempted to think of clusters as blob-like, 
well-seperated and homogeneous structures in configuration 
phase like, e. g., nano-clusters formed by agglomeration of 
atoms. The clusters which occur in high-dimensional discrete 
solution spaces are yet of a completely different nature in that 
they are more like fragmented and interweaved structures with 
lots of dead-ends, loops and holes which makes it difficult to 
speak of spatially seperated clusters. 

The existence of a clustered phase has been proven rigor- 
ously for K > 8 111 811 . In the language of statistical physics this 
clustering corresponds to one-step replica-symmetry breaking 
(1-RSB) 1I20I, |2 ill . A further substructure in terms of another 
clustering of solutions taken from one cluster, giving a hierar- 
chical structure of clusters is suspected where 1-RSB becomes 
instable and higher steps of replica symmetry breaking occur 
JH. 

What makes cluster phenomena interesting from the algo- 
rithmic point of view is the question if (and if so in what way) 
clustering has an influence on the performance of local search 
heuristics. Usually it is assumed, that the existence of many 
clusters is an indication for a complicated "rugged" energy 
landscape, which then also gives rise to many local minima, 
hindering the performance of local search heuristics II22I1 . In 
the same way, but with a slightly different focus, Krzakala et 
al. d 1 711 propose that the appearance of locally frozen variables 
in clusters is responsible for the slowdown of heuristic algo- 
rithms close to the SAT-UNSAT threshold. A locally frozen 
variable is a variable which takes the same value over all so- 
lutions belonging to one cluster. A cluster containing at least 
one frozen variable is called frozen. One defines the freez- 
ing transition <Xf as the smallest value of a above which all 
solutions belong to frozen clusters. 

To clarify the influence of phase transitions on the aver- 
age computational hardness one can study the performance of 
stochastic algorithms as a function of the control parameter 
a. Of particular interest is the algorithm-dependent value of a 
up to which an algorithm shows linear-time performance, and 
compare this to threshold values of a 12311 . Studies of stochas- 
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tic algorithms such as ASAT JH, WalkSat JH and ChainSat 
rf25tl have shown however that those algorithms have linear be- 
haviour up to values considerably beyond the clustering tran- 
sition. This suggests that the cluster transition has no impact 
on the performance of local search algorithms as long as there 
are precautions against entropic traps. It is remarkable that 
ChainSat has this behaviour although it is greedy "in a weak 
sense" as it never allows steps which increase the number of 
unsatisfied clauses. Naively one would therefore expect it to 
get trapped very easily in local minima. The authors of 02511 
interpret this as evidence for the belief that true local minima 
are very rare in high-dimensional search spaces. These results 
could also indicate that indeed it is more the (non-) existence 
of frozen clusters, which is responsible for the performance of 
local-search algorithms. 

Nevertheless, for small instances there are always some 
frozen variables. Therefore in l26ll a different notion of frozen 
clusters via the whitening core is used. There one looks, for 
each solution, interatively for variables which can be flipped 
since they appear only in clauses satisfied by other variables 
or which contain variables already detected in the whitening 
core. The position of this freezing transition was then cal- 
culated by exact enumeration and clustering of all solutions 
for sufficiently small system sizes and is expected to he at 
a,f ■ = 4.254 close to but below the satisfiability transition. 
In a similar way, rt27ll finds a cluster condensation transition 
in the solutions generated by ASAT very close to a s , again 
these results rely on a non-uniform sampling of the solutions. 
Anyway, these results are compatible with the observation of 
a good performance of local-search algorithms close to the 
threshold a s . 



C. Algorithmic treatment of SAT 

Algorithms for SAT include a broad spectrum, both 
stochastic and exact, from simple and straight-forward algo- 
rithms like RandomWalksat lEfland WalkSat to com- 
plex algorithms like DPLL IBlll and message-passing algo- 
rithms such as Belief Propagation (BP) and Survey Propaga- 
tion (SP) 02lll . For small systems exact enumeration of all 
solutions is possible using one of the numerous standard al- 
gorithms i32l l33ll such as the aforementioned DPLL. It can 
be shown IB4I1 that deterministic algorithms have longest av- 
erage run times close to oc v reflecting the difficulty of deciding 
whether a given SAT formula is satisfiable or not. The prob- 
lem with exact enumeration is that it is limited to small sys- 
tems due to hardware restrictions, especially because of the 
memory needed to store the huge number of configurations, 
as the number of solutions grows exponentially with system 
size. Furthermore the number of solutions is not a continu- 
ous function when crossing the satisfiability threshold, but it 
drops from a finite number to zero. This corresponds to a non- 
zero entropy at the phase transition, the entropy per variable 
grows approximately linearly with decreasing a 1351 l36ll . In 
turn this means that even very close to the satisfiability thresh- 
old the number of solution grows exponentially, and quickly 
becomes so large that it is not feasible to enumerate all solu- 



tions even in this regime. From counting all solutions using 
DPLL for systems up to N = 160 we can estimate the solu- 
tion entropy s = S/N near the phase boundary a = 4.25 to 
be roughly between 0.11 and 0.12 which gives more than 10 8 
solutions already for N > 160, thus taking up at least 2 GB. 

To overcome these limitations one can turn to stochastic 
algorithms which, starting at an arbitrary configuration, do 
successive changes either completely randomly or based on 
a heuristic evaluating information about the local configura- 
tional neighbourhood. Stochastic algorithms are not guaran- 
teed to find a solution, even if solutions do exist, but they can 
be significantly faster than deterministic algorithms. It is thus 
possible to obtain solutions for much larger systems, but on 
the other hand stochastic algorithms can never prove that there 
is no solution, tests for unsolvability can only be done by us- 
ing deterministic algorithms. 

In this paper we study the cluster structure numerically 
for K = 3, which requires unbiased sampling of the solution 
space. Different types of sampling algorithms are studied and 
shown to be biased. We therefore present an algorithm which 
uses a different approach to create a survey of the cluster struc- 
ture of Satisfiability instances from which it is then possible 
to derive unbiased samples. It is an improvement on the "Bal- 
listic Search" algorithm which has originally been applied to 
spin glasses II37I439I1 . The main advantage of the algorithm is 
that it is able to provide an overview of the cluster structure of 
the solution space without having to enumerate all solutions 
which is no longer possible already for moderate numbers of 
variables. 



III. SAMPLING ALGORITHMS 
A. Bias in SLS algorithms 

If stochastic local search (SLS) algorithms found all solu- 
tions with the same probability one could use them directly to 
probe the solution space. Unfortunately, this is not the case as 
we will see in this section. Later on in this article we will use 
ASAT as solution generator so we use it here for an exemplary 
presentation of the bias in SLS algorithms. 

ASAT is a simplified variant of Focused Metropolis Search 
1I24I1 and was first described in 2006 i23ll . It starts at a ran- 
dom configuration and in each step picks a variable from an 
unsatisfied clause. This variable is flipped if either this de- 
creases the number of unsatisfied clauses, or otherwise with 
a constant probability p which is a tuning parameter of the 
algorithm. ASAT has run times linear in the system size at 
least up to a = 4.21 on 3-SAT. For instances of moderate size 
like those we study here, it can well be used beyond this point 
& 

The test procedure is very simple: For a randomly cho- 
sen small instance we run ASAT again and again starting 
each time from a different randomly chosen configuration and 
count how often each solution returned by ASAT is found. If 
there were no bias we would expect the histogram of solu- 
tion multiplicities to be flat except for statistical fluctuations 
around a plateau value, i. e. the histogram should resemble a 
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solutions ordered by multiplicity 



FIG. 1: Multiplicities of solutions found by AS AT in 10 6 runs for 
a randomly chosen instance with N = 50, a = 4.0, compared to an 
unbiased distribution. Inset: Multiplicities of the ASAT solutions 
after an additional T = 0-MC step with 10 sweeps. 



Gaussian error function. 

Fig. Q] shows the resulting histogram, in comparison to a 
histogram filled with the same number of random integers 
drawn from a truely flat distribution over the range corre- 
sponding to the number of solutions of the SAT instance show- 
ing what the distribution should look like if there were no bias. 
Clearly there is a strong bias favouring some solutions over 
others. To quantify the deviations we use a % 2 test, and calcu- 
late the p-value giving the probability that an unbiased sam- 
pling process yielded a sample deviating at least as much than 
the one at hand. The p-values numerically are smaller than 
jq-323 (J e ^ t j le resolution of our double numbers). 

To test whether this bias can be corrected in a simple way, 
we did a further check, where instead of using the solutions 
returned by ASAT directly, for each solution found by ASAT 
a solution from the same cluster was generated using a T = 
Monte Carlo (MC) search starting at the ASAT result. The 
outcome of this modification is shown in the inset of Fig.[T]for 
the same SAT instance as before. The distribution now clearly 
has 5 plateaus corresponding to the 5 clusters of the solution 
space and looks much flatter but exhibits still some bias. One 
sees that most of the ASAT solutions stem from the smallest 
cluster, hence the sampling does not respect the cluster size. 
Hence, ths bias can be strongly decreased by additional T = 
MC simulations, but not completely. Further checks showed 
that the bias persists independently of the system size. 

Since we want to study clustering properties of the solution 
ensemble we need to remove the bias completely and sample 
solutions in proportion to the cluster sizes. To ensure this, we 
will perform reweighting using the Ballistic-Search algorithm 
as described in section [IV] Before we come to the Ballis- 
tic Search, we will show in the next section that MCMCMC, 
another important sampling method, fails on sampling SAT 
solutions uniformly as well. 



B. Bias in MCMCMC 

The Metropolis-Coupled Markov Chain Monte Carlo 
(MCMCMC) method, first proposed in 1991 by Geyer |40ll . 
also known as Parallel Tempering [41, 42], is a powerful and 
versatile tool, commonly used in biophysics and statistical 
physics to perform equilibrium simulations and to generate 
unbiased samples in large configuration spaces. MCMCMC 
uses a set of replicas of single instances, simulated in par- 
allel at different temperatures and linked by global updates 
in which replicas are swapped pair-wise with an acceptance 
probability depending on their energy difference and temper- 
ature spacing (Metropolis-Hastings criterion), thus facilitating 
the tunneling through barriers seperating local minima of the 
phase space H43I1 . 

To study the performance of MCMCMC on SAT we em- 
ploy a histogram test similar to the one described in sec- 
tion IIII Al for the performance of ASAT. For several values 
of a = 1.00... 4. 25 scattered over the satisfiable phase, the 
number of variables N is chosen such that expected number 
of solutions is 1000. This, e. g., results for the smallest value 
of a considered here in a system size N = 14, while for the 
highest value of a, N = 50 is feasible. 

We apply a straight-forward implementation of MCMCMC 
to a set of 50 instances for each value of the control parame- 
ter a, where we use 15 temperatures, the lowest, at which the 
samples are taken, being initially 7b =0.1, the highest such 
that the corresponding energy is found to be approximately 
M2~ K which is the expected energy of a completely random 
configuration. Every 1000 steps the temperatures are adjusted 
to drive the replica exchange rate between neighboring tem- 
peratures towards 50 lowest temperature fixed at 0.1. The 
procedure chosen to adjust the temperatures leads to a distri- 
bution of temperatures where for the lowest temperatures the 
exchange rates indeed reach 50 whereas the highest tempera- 
tures all gather in the random phase. This can be seen as an 
indication that the number of temperatures used is sufficient 
to allow the replicas to travel between the constrains, i. e. the 
highest temperatures are indeed located in the "paramagnetic 
phase". We take one sample every second sweep to generate 
a total of 10 6 samples. Only successful sampling steps are 
counted, i. e. those where the energy of the configuration at 
7q is zero. 

The histograms with the resulting distribution look pretty 
close to those drawn from a flat distribution (not shown here). 
We again use the p-values obtained from a % 2 test to quantify 
the deviations and find that in most cases MCMCMC gives 
reasonably flat distributions, hence this method appears to ex- 
hibit on the first sight a much lower sampling bias. Never- 
theless, there also are a number of histograms having a sig- 
nificant bias exhibiting very small /^-values. The higher the 
value of the control parameter a is chosen, the larger becomes 
the spread of the distribution of p-values towards extremely 
small values. This can be seen from Fig. [2] where the distri- 
bution of p-values, integrated over all system sizes and values 
of the control parameter, is shown for the two cases where the 
instances exhibit one or more than one cluster of solutions, 
respectively. (The number of clusters can easily be calculated 



5 



0.2 



3 



0.4 0.6 
p value 



FIG. 2: Bias of MCMCMC: Dependence of p-value on the number 
of clusters, integrated over all N and a. Dashed line: p-values for 
instances with only one cluster, solid line: p-values for instances 
with more than one cluster. 



exactly for these small instances.) In the case that the solution 
space consists of only one connected component, the p-valne 
distribution is flat showing that MCMCMC works unbiased as 
expected. The presence of clustering on the other hand leads 
to a bias or imbalance in the sampling process resulting in a 
strong peak at small p-values. 

Since we are interested in particular in those instances 
which exhibit many clusters, MCMCMC turns out to be not 
suitable as well, since all instances have to be sampled cor- 
rectly. Note that for larger system sizes, the number of in- 
stances having just one cluster, where MCMCMC seems to 
work well, will strongly decrease. Hence, for large system 
sizes, MCMCMC will exhibit a bias for basically all instances 
of interest. To create an unbiased sample we need a different 
method which will be presented in section HVl 



IV. BALLISTIC SEARCH 



Here, as mentioned in Sec. Ill Bl we are using the neighbor- 
based definition of clusters: two solutions are considered to 
be in the same clusters, if there exists a path in solution space 
consisting of single-variable flips. 

We use Ballistic Search, which has been introduced in the 
year 2000 as a method for studying ground-state properties of 
spin glasses [37]. The approach is able to provide a survey of 
the cluster landscape using stochastic algorithms, in particular 
without the need to enumerate all ground states as it is usually 
necessary when one aims at clustering. The sheer number of 
ground states forbids exact enumeration when studying spin 
glasses, and, as mentioned above, the same holds for SAT. We 
therefore use this method which relies on generating a survey 
of the most important clusters. 

The survey consists of a set A = {A,} where each ele- 
ment A, = ({c^ 1 ^},£W) represents one cluster and consists of 




FIG. 3: Between those 6 solutions (black circles) from the same clus- 
ter (light grey) Ballistic Search has found only 2 paths. The apparent 
number of clusters is 4. We need to increase the density of solutions 
to make Ballistic Search more efficient. 



a (small) set of solutions {cj } from the cluster i and an esti- 
mate £W of the size of cluster i. The survey should cover all 
clusters, or at least all but those which are negligibly small. 
One can then sample the whole solution space with correct 
weights by generating the desired number of solution samples 
from the representative sets of solutions for each cluster ac- 
cording to the respective cluster sizes. 

Two main ingredients form the basis of the Ballistic Search 
algorithm. Firstly, the above described data structure storing 
small sets of representative solutions for each cluster instead 
all solutions. Secondly, a "ballistic path search" is used to 
analyse the cluster space and generate the survey from a given 
set of solutions. The basic operation of this procedure is that 
we have to determine for any given pair c a , ct, of solutions, 
whether they belong to the same cluster. This has to work un- 
der the assumption that for the case of c a , c\, belonging to the 
same cluster, a complete nearest-neighbor path of solutions 
between c a and c\, is not contained in the set of already found 
solutions. Instead, one searches stochastically for paths be- 
tween c a and Ch by starting at one solution and subsequently 
changing a randomly chosen free variable. A variable is called 
free if its value can be changed without violating any con- 
straint, so that one never leaves the solution cluster. This is 
repeated until either the target solution is reached or no free 
variable is left, because as an additional stop condition every 
variable shall be touched at most once. (Because of this ad- 
ditional constraint the path search is called "ballistic".) This 
implies that in a successful ballistic path search the number 
of steps taken is always equal to the Hamming distance of the 
solutions, i. e., the number of variables in which the two solu- 
tions differ. 

Figures|3]and|4]give a graphical description of how the Bal- 
listic Search algorithm works. We start with some randomly 
generated solutions depicted as black circles in Fig. [3] all 
of which belong to the same cluster which is drawn in grey 
in a 2-dimensional cartoon of the N-dimensional configura- 
tion space. For the sake of simplicity we assume here that 
all solutions belong to the same cluster, the generalization to 
more than one cluster is obvious. Running the ballistic path 
search we find that some of the solutions can be connected by 
paths drawn as lines in the picture, i. e., for these solutions the 
algorithms correctly finds that they belong to the same clus- 
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FIG. 4: Adding solutions (grey) has yielded a correct identification 
of the cluster, because more paths (grey) have been found, now con- 
necting all solutions. 




FIG. 5: This cluster has a more complex structure than the one in 
Fig. [4] illustrated by the additional holes. Here adding even more 
new solutions (grey) than before does not work. Still not all solutions 
are recognized as belonging to the same cluster. 



ter. The problem is that for low solution densities the average 
distance between solutions is large, and the efficiency of the 
ballistic -path search strongly decreases with larger distances 
l37tl . Therefore we only find a few paths and the apparent 
number of clusters in this example is larger than its true value, 
cf. Fig. [3] What we need to do is to increase the number 
of solutions by rerunning ASAT. For few added solutions, the 
measured number of clusters will increase, since only few ad- 
ditional paths within clusters are detected, less than the num- 
ber of added solutions. When generating even more solutions, 
we will find that the apparent number of clusters at some point 
no longer increases, but instead it decreases as more and more 
paths between solutions are found, until finally all solutions 
are correctly assigned to the same cluster as shown in Fig. [4] 



A. Ballistic networking 

Our studies have shown that the simple ballistic-path search 
algorithm has very low efficiency when applied to Satisfiabil- 
ity which can be attributed to a high complexity of the solu- 
tions space or large sizes of the clusters. We therefore devel- 
oped a refinement of the algorithm named "ballistic network- 
ing", which is a very general extension of the ballistic path 
search so that it can readily be applied to other problems. 

The idea of the algorithmic refinement is to increase the 
probability of identifying two solutions, origin c a and target 



FIG. 6: Ballistic Networking improves on the result of Ballistic 
Search by not adding solutions randomly, but adding solutions from 
the same cluster using al = MC search (arrows). Now all solu- 
tions are found to belong to one cluster. 



Cb, as belonging to the same cluster using ballistic path search 
by again increasing the number of solutions. Instead of using 
ASAT to generate more solutions we generate In^M additional 
solutions by performing independent T = MC simulations 
starting at c„ and q,, respectively. Hence, we are sure that the 
additional solutions belong to the same cluster as their respec- 
tive "parent" solution. We then try to find connections using 
the ordinary ballistic path search between all (n a dd + l) 2 pairs 
of solutions, where one solution belongs to the origin and one 
to the target. If at least one path is found, it is clear that c a and 
ct belong to the same cluster. We apply this test to all pairs 
of solutions which have not yet been found to belong to the 
same cluster. An artist's view of this improvement is given in 
figures[5]and[6] Fig. |5]shows how the standard ballistic search 
fails due to a more complex structure of the cluster, although 
even more solutions than before have been used. In Fig. |6]the 
solutions found by the T — MC search are drawn as points 
connected to their parent solution by arrows. We can see that 
the number of successful ballistic path searches (gray lines) 
does not have to be very high, but still is enough to correctly 
identify the cluster. 

Indeed this procedure improves the performance of the 
search so much that it outweighs the additional effort of hav- 
ing to carry out (n a dd + I) 2 ballistic path searches instead of 
one. Fig. |7]shows a comparison of the performance of ballis- 
tic path search without and with additional solutions. The case 
"+0" corresponds to the bare ballistic path search. The hori- 
zontal axis shows the number of ballistic path searches which 
have to be carried out in the worst case, i. e., when no con- 
necting path is found. We found 5 < n a dd < 10 to be a suitable 
range for the system sizes under study. 

When creating the additional solutions to test whether a pair 
of solutions c a , c\, belongs to the same cluster, to improve the 
success probability, one can think of introducing a bias into 
the T = MC search which pushes the additional solutions 
derived from the first solution c a closer to the second solution 
Cb and vice versa. Indeed we found that such a bias has a pos- 
itive influence on the success probability of the ballistic path 
search. Yet we did not use this bias in the implementation, 
because the positive influence comes at a high cost. For each 
pair of solutions to be tested a dedicated biased set of addi- 
tional solutions has to be generated which cannot be reused 
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FIG. 7: Comparison of ballistic path search without and with addi- 
tional solutions ("ballistic networking") for N = 128 and a = 3.0. 
We show the probability for finding a path between two solutions 
generated from the same cluster, as function of the total number of 
ballistic path searches between all pairs of parent and children con- 
figurations averaged over 1000 runs. The case "+0" corresponds to 
the original ballistic path search. 




actual cluster size 

FIG. 8: Comparison of the actual number of configurations per clus- 
ter to the number estimated by Monte Carlo integration. Each point 
corresponds to one cluster. 



when comparing either c a or q, to a third solution. The neces- 
sary computational effort for generating each time new biased 
configurations by far outweighs the positive effect of the bias. 

The second part of the cluster survey A consists of the sizes 
£W of the clusters. An exact calculation of the cluster size is 
possible, but takes too long, since it typically grows exponen- 
tially with N. We have therefore examined several different 
estimation methods with respect to their reliability in giving 
a correct estimate for the cluster size by comparing the esti- 
mated cluster size to the exact cluster size on a random ensem- 
ble of clusters for different values of a and small system sizes 
N. The best method known to us has been found to be the es- 
timation of the cluster size using a Monte Carlo integration as 
it has been used in 14411 in an application to spin glasses. Fig. 



[8] shows a comparison of the actual number of configurations 
in one cluster to the number as estimated by Monte Carlo in- 
tegration. We used several combinations of N and a where 
the total number of solutions (over all clusters) did not exceed 
5 • 10 6 , such that all clusters could be calculated exactly, and 
afterwards for each cluster the MC estimation was run. 



B. Implementation 

Combining ballistic networking and the cluster size estima- 
tion the full algorithm is comprised of two alternating steps. 
The first step is to generate a given number (of the order of 
1000) of solutions of the Satisfiability instance at hand using 
the AS AT algorithm. The tuning parameter of AS AT is chosen 
to be p = 0.21 which is the optimal value as given in lr23"ll . In 
the second step Ballistic Networking of the solutions found by 
ASAT is done as described above to create the cluster survey, 
and then the sizes of the clusters are estimated. Afterwards 
ASAT is run again and another set of new solutions is created. 
The cluster survey is then updated using Ballistic Networking 
on the new solutions and the solutions representing the so far 
found clusters in the existing survey. Here new clusters may 
be found, and if so, their size is estimated. This is repeated 
until the cluster survey is considered complete, i. e., no more 
relevant clusters are found. 

From the cluster survey for each instance a set it of un- 
biased solutions can be generated using the cluster-size esti- 
mates. For each solution to be generated for a given instance, 
first a cluster from the survey is selected with a probability 
proportional to the cluster size. One solution is selected from 
the set of representative solutions and starting from this solu- 
tion a T = MC search is performed finally giving the solu- 
tion to be used in the analysis. 

Defining a good stopping criterion is a crucial point of the 
algorithm. As the cluster number in Satisfiability can be rather 
large, we decided not to generate all clusters, but all except 
for those which contain only a neglegibly number of solu- 
tions. For this purpose we monitor the total cluster weight 
We run the algorithm until the total cluster weight has 
not increased by more than 0.5 % over the last half of solu- 
tions included in the clustering process. We store the order in 
which the solutions have been generated by ASAT and label 
each cluster with a number telling the position of the earliest 
solution which has been found to belong to this cluster. 

When trying to optimize the number of new solutions added 
in each round one has to consider two competing effects: 
On the one hand adding solutions — as in ordinary Ballistic 
Search — may reveal that two clusters actually are parts of the 
same cluster, connected maybe by only a small path in config- 
uration space which has been too hard to find with fewer so- 
lutions. On the other hand increasing the number of solutions 
makes Ballistic Networking slower and, even worse, increases 
the probability of false new clusters which in turn can lead to 
an on-going increase of the cluster number and total cluster 
weight and thus to a failure of the stopping criterion. 

The system sizes which can be reached using the method 
described above depend, of course, on the control parameter 
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FIG. 9: Dependence of p pat h of ballistic path search on a and djjam 
for N = 128. As the typical distance of two solutions depends on a, 
so do the ranges shown in the plot. For small distances dham < 15, 
a path is usually found, while for large distances, the probability de- 
pends on the value of the control parameter a, which can be under- 
stood by a more and more complex structure of the solution clusters. 



a. For small a all solutions are contained in only one large 
cluster where there are many possible paths between config- 
urations so that Ballistic Search is very efficient and system 
sizes of a few hundred variables are possible. For high a val- 
ues in the solvable phase, the number of solutions is small so 
that in this regime Ballistic Search still is rather efficient due 
to the small extent of the clusters and relatively large system 
can be done. 

Fig. [9] shows the dependence of the success probability of 
ballistic path search on the Hamming distance dham between 
the configurations for N = 128, for different values of the con- 
trol parameter a. Up to a — 3.8 the probability decreases 
strongly with increasing a, as the clusters develop more holes. 
Above this point the curve is approximately independent of a. 
We also find that the probability decreases weakly with in- 
creasing system size (not shown). The fact that the average 
distance between solutions decreases with a makes Ballistic 
Networking most difficult in the intermediate regime around 
a « 3.3. Here the number of ballistic path searchs needed to 
find a connection between two solutions from the same clus- 
ter is highest. The cluster structure seems to be such that there 
are many "dead ends" in which the search may get stuck. To- 
gether with the high number of clusters, which enters quadrat- 
ically in the running time, this limits the reachable system 
size. All in all, Satisfiability instances of up to N = 144 vari- 
ables were doable in reasonable time over the whole range of 
interest 3.0 < a < 4.2 while for smaller intervals of the control 
parameter, we also studied N = 256. 



V. RESULTS 

Here, we study the behavior of random 3-SAT instances as 
a function of the parameter a. This is meant in the sense that 
we generate an instance using a given number N of variables 



and a set of (arbitrarily ordered) clauses C m (m = I , . . . ,M max ). 
We chose M max = 0C max A, where 0C max is the largest value of 
the control parameter we want to consider. We can study the 
behavior of each instance as function of M < M max by con- 
sidering each time exactly the clauses C m for m = 1,... ,M. 
Also, we can average over these distances for each value of 
the control parameter. 



A. Hierarchical clustering 

For the analysis |45] of the behavior of 3-SAT as a function 
of the control parameter a, we start by looking at the hierar- 
chical structure of a set il of solutions sampled for a typica l 
3-SAT instance. We have used "Ward's algorithm" j46l l47ll . 
an agglomerative hierarchical matrix updating algorithm, on 
the set if to extract a hierarchical clustering from which we 
can then draw a visual representation of the solution space. 

Ward's algorithm has been applied in many different fields 
ranging from RNA secondary structures over optimization 
problems to spin glasses 10, [HI HH HH] . It is an iterative pro- 
cedure where initially each configuration comprises a single 
item cluster. In each step those two clusters are merged which 
have minimal distance with respect to an effective distance 
measure chosen such that the sum of the variances in each 
cluster is minimized. After each merger, the distances of the 
remaining clusters to the new cluster have to be calculated, 
for details see, e. g., Ref. [47]. Finally, one reorders the con- 
figurations according to the hierarchy obtained in the iterative 
merging process, and draws a color-coded visualization of the 
distance matrix. 

Next, we present some results for a typical instance. We 
chose one which exhibits its SAT-UNSAT transition close to 
the numerical estimate of the ensemble average a s = 4.267 
given in II 1611 . Fig. [10] shows the color-coded distance matri- 
ces and the dendrogram which were generated for three dif- 
ferent values of a . The difference in the solution landscape 
and cluster structure between the phases is clearly visible. For 
low a the Ward matrix is featureless and homogeneously grey. 
All solutions belong to one single cluster and the phase space 
shows no specific features. In the intermediate range one sees 
box-like structures along the diagonal in a darker grey. These 
correspond to clusters, because darker means smaller Ham- 
ming distance and the solutions inside a cluster are closer to 
each other than to other solutions. Some of these boxes show 
a substructure which can be interpreted as the solutions from 
this cluster themselves forming sub-clusters. This is consis- 
tent with the theoretical prediction of replica-symmetry break- 
ing beyond 1-RSB in the intermediate a range 112211 . Never- 
theless, as mentioned in the introduction, it is to be expected 
that most of the clusters are not relevant in the thermodynamic 
limit and a small number of clusters contains almost all solu- 
tions. For higher values of a the substructures inside the clus- 
ters become washed out whereas the first-level cluster struc- 
ture becomes more pronounced as the cluster become smaller. 
In the replica symmetry breaking framework this would be in- 
terpreted as a vanishing of higher level RSB above a certain 
threshold, but this cannot be deduced from looking at single 
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FIG. 10: The hierarchical structure resulting from Ward's algorithm 
visualized both as tree structure (dendrogram) and distance matrix, 
for N = 256 and a = 1 .00 (top), 4.00 and 4.25 (bottom). Darker grey 
scales correspond to smaller distances. 



instances, of course. 



FIG. 1 1 : Complexity c as a function of a and for several system sizes 
between 32 and 144. Lines have been drawn to guide the eye. The 
error bars give statistical errors. 



B. Averaged quantities 

The complexity c — i logA^ £ is defined as the logarithm of 
the number of clusters normalized to the system size. Fig. [TT] 
shows the complexity as a function of a and for system sizes 
up to = 256 variables averaged over 200-500 instances for 
each value of a and each value of N. The number of clusters 
was taken directly from the cluster surveys created using the 
Ballistic Networking method described above. 

For the "easy" part of the satisfiable phase, where the value 
of the control parameter a is small, there is only one clus- 
ter, thus the complexity is zero. In an intermediate range the 
number of clusters grows peaking at a value which is strongly 
affected by finite-size effects and then becomes smaller again. 
This behaviour reflects the theoretical prediction of one single 
cluster in the low a regime "crumbling" into smaller pieces 
when a is increased and the clustered phase is reached. For 
even higher a the vanishing of solutions leads to the disap- 
pearance of clusters and the cluster number decreases again. 
Clearly, the peak of the complexity curves seems to converge 
towards a c with increasing system size, in accordance with 
the analytic predictions IT17I1 . 

Looking at this plot one has to keep in mind that the stop- 
ping criterion used in the algorithm is based on the number 
of solutions covered by the clusters that were found so far. In 
a phase with a large number of small clusters we will miss 
small clusters if they only comprise a negligible part of the 
solutions (in the sense of the stopping criterion) and therefore 
underestimate the number of clusters. It is therefore natural 
that the complexity found here is lower than the one given in 
ll26ll . After all the complexity shown in the graph is only a 
lower bound for the true complexity respecting all clusters. In 
phases dominated by few and large clusters its value should 
nevertheless be close to the true value. 

Fig- Q~2] displays the fraction of the solutions contained in 
the largest cluster. For a < a c this value seems to increase 
with growing system size. At a = tt c it exhibits a minimum, 
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FIG. 12: Fraction of the weight of the largest cluster with respect to 
the total weight of all clusters. 



while for a > a £ it decreases slightly with growing system 
size, but it is larger than the value found right at a £ . These 
results are also compatible with the analytical prediction lfl7ll . 
which states that only for a range a > a c more than one cluster 
is relevant in the thermodynamic limit. Nevertheless, we can- 
not deduce from the data, since system sizes we can reach are 
rather limited, whether for all values of a < a c this growing 
fraction converges to one or to a smaller value. 

Next, we have a closer look at the average structure of the 
solution space. As mentioned in the discussion of the Ward 
matrices, solutions belonging to the same cluster are more 
similar to each other, i. e., closer in terms of the Hamming dis- 
tance, than pairs of solutions which belong to different clus- 
ters. The cluster structure is thus reflected in the set of all 
pairwise overlaps, where the overlap r,y of two solutions i and 

j for which the Boolean variables take the values given by x$ 

and xjP is defined as rij := Yl^=\ 8(*n ,x}P)lN. Fig. [T3l shows 
the overlap distribution for a = 4.0 and several values of TV. 
For each system size at least 1000 instances have been pro- 
cessed with the algorithm described in section IIVBI and 500 
solutions have been generated from each cluster survey. 

Two peaks are visible. One peak is lying close to (r) = 1 
and due to the overlap of solutions belonging to the same clus- 
ter. With larger system size it moves slightly to lower (r) val- 
ues and becomes sharper. The second peak at about (r) = 0.7 
is not discernible for the smallest system size, but only evolv- 
ing with larger system sizes and only visible weakly against an 
also growing background. Note that a pure two-peak structure 
would correspond to the picture of one-step broken replica 
symmetry (1-RSB) ll49ll . Nevertheless, the result is not fully 
clear here, since in addition to the peaks, there is also a con- 
tious part between the two peaks. On the other hand, it is clear 
that the overlap converges to zero for values of the overlap 
smaller than 0.5. This speaks against a full-level RSB struc- 
ture of the solutions space. Note that we have found similar 
results for other values of the control parameter a c < a < U s 
(not shown). 
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FIG. 13: Overlap (r) of solutions in clusters found by Ballistic Net- 
working, for a = 4.0. For (r) < 0.4 the curves are essentially zero. 



C. Freezing transition 

To complete the picture we also studied the freezing tran- 
sition, which as mentioned in section HI Bl is defined as the 
smallest a above which all solutions belong to frozen clusters 
and has been found to lie at a/ = 4.254. 

To check directly whether a cluster contains frozen vari- 
ables, we need to generate and compare all solutions from 
this cluster, therefore cluster surveys do not help here. Us- 
ing exact algorithms we find that for the system sizes we can 
reach, for all a near the SAT-UNSAT transition there are al- 
ways frozen variables in all clusters. This is probably due to 
too small systems sizes. 

Thus, we followed a different approach. For each in- 
stance, taken at a = 4.20 and 4.25 and for system sizes up 
to N = 2048, we generated a solution using ASAT, which be- 
longs with high probability to the largest cluster. Then we 
performed a very long T = simulation starting from this so- 
lution, and measured the fraction pf 10Z en of variables which 
have never flipped while performing this random walk inside 
the solution cluster. We extrapolated this fraction to a large 
number of MC steps, yielding Pfy ozen , see Fig. [14] With in- 
creasing system size, Pfy ozsn seems to converge to zero, see 
inset of Fig. [14] Hence for a = 4.20 and a = 4.25 the largest 
clusters seems to contain no frozen variables in the thermo- 
dynamic limit. This is compatible with a.f = 4.254, meaning 
that in the thermodynamic limit no frozen clusters occur be- 
low this value of a. 



VI. CONCLUSION 

In this work we have shown that stochastic local search al- 
gorithms cannot be expected to produce correctly weighted 
samples of the solution space of Satisfiability. The same holds 
true for MCMCMC which is widely used to sample config- 
uration spaces in many fields of application, when the SAT 
solutions are spread over several clusters. 
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FIG. 14: Fraction pf mmn (t) of never-flipped variables for the 
"largest" (i. e. first detected) cluster during T = MC simulations for 
a = 4.20 and different system sizes N. The solid lines indicate fits to 
functions of the form pf rozen (t) = Pf mzea +at~P Inset: Extrapolated 
values Pff 0zen as function of system size N in a double-logarithmic 
scale, for a = 4.20 and a = 4.25. The "?" marks a point where the 
extrapolation failed and an upper limit was estimated by eye. 



A new type of algorithm has been presented and used for 
studying clustering phenomena in the solution landscape of 
the Satisfiability problem. It is an improved version of the 
Ballistic Search algorithm which has been successfully used 
for studying spin-glasses. Its guiding principle is to gener- 
ate a survey of clusters of solutions represented by small sets 
of solutions rather than enumerating and clustering all solu- 
tions which is unfeasible already for moderate system sizes. 
By using a different approach, Ballistic Networking, in the 
reconstructing process of the cluster structure the efficiency 
of the Ballistic Search could be improved so that its perfor- 
mance becomes reasonably high when used on Satisfiability. 
The method presented here is general enough to be suitable for 
many other problems. Of course, it would be natural to study 
Satisfiability for K > 3 using Ballistic Networking, but the ef- 



ficiency for ballistic path search seems to be still much lower 
than for the case of K = 3 which sets very restricting limits 
on the system sizes which can be reached. Nevertheless, the 
approach presented here should be useful for many disordered 
systems like other types of combinatorial optimization prob- 
lems. 

In the case of Satisfiability, the range of low values of a 
(where many solutions exist but belong to only one cluster) 
can be studied by MCMCMC. Furthermore, the case of high 
values of a close to the SAT-UNSAT transition (fewer solu- 
tions contained in several clusters) can be studied using exact 
enumeration of all solutions. In contrast, the algorithm pre- 
sented here allows to study the full satisfiable phase, but it is 
limited to moderate system sizes in the intermediate a range. 
Nevertheless it is the only reliable method to generate unbi- 
ased samples in this regime. 

Using the method described here the ensemble properties 
of Satisfiability with moderate system size could be studied 
and analytic predictions about the cluster structure could be 
tested. To this aim we first did a visual inspection of the 
cluster landscape using a graphical representation in terms of 
Ward distance matrices. These show the expected structural 
differences of the different phases of Satisfiability. Further- 
more we had a look at the complexity measure over the whole 
a spectrum in the easy phase, and the overlap distribution of 
the solutions for particular values of a. Our findings are in 
good agreement with the theoretical predictions and previous 
numerical studies using other methods. 
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